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We propose an efficient Monte Carlo algorithm for the off-lattice simulation of dense hard sphere polymer 
melts using cluster moves, called event chains, which allow for a reject ion-free treatment of the excluded 
volume. Event chains also allow for an efficient preparation of initial configurations in polymer melts. We 
parallelize the event chain Monte Carlo algorithm to further increase simulation speeds and suggest additional 
local topology-changing moves (“swap” moves) to accelerate equilibration. By comparison with other Monte 
Carlo and molecular dynamics simulations, we verify that the event chain algorithm reproduces the correct 
equilibrium behavior of polymer chains in the melt. By comparing intrapolymer diffusion time scales, we show 
that event chain Monte Carlo algorithms can achieve simulation speeds comparable to optimized molecular 
dynamics simulations. The event chain Monte Carlo algorithm exhibits Rouse dynamics on short time scales. 
In the absence of swap moves, we find reptation dynamics on intermediate time scales for long chains. 


I. INTRODUCTION 

Polymer melts or polymer liquids are concentrated so¬ 
lutions of long chain molecules above their glass or crys¬ 
tallization temperature. In a dense polymer melt long- 
range excluded volume interactions become screened and 
an individual polymer shows ideal behavior.® Polymer 
melts exhibit a characteristic and complex dynamical 
and rheological behavior because of entanglement effects, 
which impede chain diffusion and give rise to reptation 
dynamics of polymer chains.®® The melt state is also 
most relevant for processing and manufacturing polymer 
materials.® 

In this paper we introduce a novel Monte-Carlo (MC) 
algorithm for the off-lattice simulation of a melt of flexi¬ 
ble hard sphere p(^mer chains, which are connected by 
springs or tethers.®®^ This event chain (EC) algorithm 
allows for a much faster equilibration as compared to MC 
algorithms based on local moves. 

The simulation of polymer melts by Molecular Dynam¬ 
ics (MD) or MC simulations is a challenging problem, in 
particular, for long chains at high density, where poly¬ 
mers in the melt exhibit slow reptation and entangle¬ 
ment dynamics.® Eor chain molecules of length N, the 
entanglement time increases oc ^ which impedes the 
equilibration of long chain molecules in a melt if only lo¬ 
cal self-avoiding displacement moves of polymer segments 
are employed as in a typical off-lattice MC simulation. 
In order to reach equilibrium by such local moves, the 
system has to go through slow reptation dynamics on 
time scales between the Rouse and entanglement time. 
In MD s imula tions, such reptation dynamics has been 
observed) ^ In MC simulations, indications of repta¬ 
tion dynamics have been observ ed in lattice models^ or 
fluctuating bond lattice models.l^^IIIl Xo our knowledge. 
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reptation dynamics has not yet been observed in an off- 
lattice MC si mulation so far, where equilibration is more 
difficult 

The dynamics of MC simulations depends on the MC 
moves that are employed. Eor local MC moves, the 
polymers obey Rouse dynamics on short time scale^^^^ 
until entanglement effects eventually give rise to the 
crossover to reptation dynamics if MC moves obey the 
self-avoidance constraint This means that the re¬ 
sulting MC dynamics can resemble the actual motion 
of coarse-grained polymers, although the MC dynam¬ 
ics is not explicitly based on a real istic mic roscopic 
dynamics.^ Local MC reptation movej^ ^ * ^^ * ^^ * (slither¬ 
ing snake moves) are used to initiate reptation dynam¬ 
ics and obtain faster equilibration of a polymer melt. 
MC simulations have the general advantage that also 
non-local or collective MC moves can be introduced, 
for ex ample, chain-topology changing double-bridging 
moveswhich speed up equilibration (such moves 
can also be combined with MD simulations to equili¬ 
brate the systenPH). Dynamic properties, however, are 
no longer realistic if such topology-changing moves are 
employed. In particular, reptation dynamics will not oc¬ 
cur if chain-topology changing moves are employed. 

If polymers in a melt are modeled as bead-spring mod¬ 
els with hard sphere beadsan additional simulation 
problem arises, in particular in MC simulations. At high 
segment or monomer densities, the mean free path of seg¬ 
ments is limited and local MC di^lacement moves are 
restricted to very small step-sizes.® 

Eor hard sphere systems, non-local cluster moves rep¬ 
resent a successful strategy to overcome the problem of 
slow MC equilibration in general by reducing rejection 
rates in the dense limit. In Ref. [22] the reject ion-free 
event chain algorithm has been proposed, which coher¬ 
ently moves large clusters of particles in the form of a 
chain, and a significant speed-up in the sampling of the 
hard sphere system has been shown. The EC algorithm 
can be generalized from athermal h ard sp here systems 
to spheres with interaction potentials) ^^ * ^^ * In Ref. [25j we 
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showed that the EC algorithm can be used for simula¬ 
tions of semiflexible bead-spring polymer systems. In 
this work, we adapt the EC algorithm for the MC simu¬ 
lation of dense polymeric melts consisting of flexible hard 
sphere polymers, verify the algorithm, and benchmark its 
performance. 

The paper is structured as follows: in Section |n| we 
present our EC based MC algorithm for hard sphere 
polymer melts. In order to further improve perfor¬ 
mance we also introduce a parallelized version of the 
EC algorithnP^ and a version employing local topology¬ 
changing “swap” moves. Eurthermore, we show that 
EC moves can also be used to efficiently generate ini¬ 
tial polymer configurations for the simulation, which are 
already representative of equilibrium configurations. In 
Section |III[ we verify our algorithm by a detailed com¬ 
parison of equilibrium structural properties, such as the 
polymer shape and the end-to-end distance distribution, 
to simulation results from other MD and MC simula¬ 
tion techniques. Naturally, these results are not novel. 
Therefore, all details of this validation are presented in 
the Appendix. Einally, in Section El we benchmark 
the performance of serial and parallelized EC algorithms 
with or without swap moves against standard local MC 
schemes and against state of the art MD simulations (us¬ 
ing the LAMMPS packag^^. We use time-dependent 
mean-square displacements (MSDs) of polymer beads to 
monitor inter- and intrapolymer diffusion and use the in¬ 
trapolymer diffusion to compare the performance of all 
algorithms in terms of a polymer relaxation time. The 
EC algorithm obeys Rouse dynamics on short time scales. 
Moreover, we show that, in the absence of topology¬ 
changing swap moves and for long polymer chains, the 
EC algorithm exhibits reptation dynamics on intermedi¬ 
ate time scales, before a crossover to chain diffusion on 
the longest time scales. We end with a conclusion and 
outlook. 


II. EVENT CHAIN ALGORITHM FOR POLYMER 
MELTS 


A very fundamental model for a self-avoiding flexible 
polymer is a bead-spring model, in which all beads in¬ 
teract via an excluded volume constraint, i.e., the poly¬ 
mers consist of hard impenetrable spheres, and the beads 
in one polymer are bonded with Hookean springs. The 
spring constant has to be sufficiently large as to enforce 
the impenetrability of polymers and avoid unphysically 
large bond stretching. In summary, we have a hard 
sphere interaction between all pairs of beads 
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with the diameter a of the hard spheres, and a harmonic 
stretching energy which, for a single polymer, can be 



FIG. 1. Simulation snapshot of a polymer melt at a vol¬ 
ume fraction r] ~ 0.54. The color of the beads discriminates 
the individual polymers. To give better insight into the poly¬ 
mer melt structure, we do not wrap polymers periodically 
although we employ periodic boundary conditions. 


written as 


k 

^bonds = 2 ^ (2) 

^ i=l 

Here, k is the spring constant, — 1 is the number 
of bonds in a polymer (containing N beads), bi is the 
length of the i-th bond, and the equilibrium length of 
the bonds coincides with the hard sphere diameter a in 
our model. In our simulations, we chose bond stiffnesses 
{ka^^/ksT = 30) such that thermal bond stretching re¬ 
mains weak with (6^) ~ l.lcr. To simulate a polymer melt 
at a given density p, we generate a system of M polymers 
in a cube of edge length L = 40cr, see Eig. We employ 
periodic boundary conditions in all directions. 

Alternatively, we also consider systems of hard sphere 
polymers bonded by tethers of maximal length 6max = 
1.4(7 rather than springs.!^ 

In systems of dense hard spheres, standard Metropolis 
MC schemes based on local moves of individual spheres 
suffer from very slow sampling, as the move length is lim¬ 
ited to roughly the mean free distance between spheres. 
This has been overcome by the introduction of suit¬ 
able cluster moves, the so-called ECsP In an earlier 
work, we extended this approach to parallel computation 
and demonstrated how the EC algorithm can be applied 
to dense polymer systems.!^ Because the EC algorithm 
moves dense regions of hard spheres or polymer beads 
coherently, it also mimics the essential features of the 
actual physical dynamics on a coarse time scale such as 
diffusion of polymer bundles.l^ If the bonded beads in a 
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FIG. 2. (a) Construction of an EC cluster in a hard sphere system. The displacement vector v is distributed to all colored 

beads, which are part of the EC cluster. The final positions are shown in lighter color, (b) EC displacements for harmonic bond 
energies if the green bead is displaced, either in direction vi without hard sphere collision or in direction V 2 where it collides. In 
both cases, the red bead becomes the next pivot bead. Eor explanations, see main text, (c) Reflection of a spring-triggered EC 
necessary because of the decomposition into parallel simulations cells. The EC cannot be transferred to the gray bead, which 
is rendered immotile. The pivot bead (green) does not change, and the propagation direction is reflected. Eor explanations, 
see main text. 


polymer interact only via pair potentials, such as in the 
present example of hard sphere bead-spring polymers, a 
completely reject ion-fr ee simulation solely based on EC 
moves is possible 

In the EC algorithm, we first choose a total displace¬ 
ment length which is the same for all EC moves. For 
hard spheres, each EC move is constructed according to 
the following rule, which is also illustrated in Fig. [^a). 

1. Select the starting pivot bead for the EC move and 
a direction (which we call v) randomly. Initially, 
the remaining displacement is Xmax = 

2. Evaluate the largest possible displacement x < 
^max of the pivot bead in the chosen direction be¬ 
fore it touches another bead. Move the pivot bead 
by X. 

3. Continue the EC move at the new pivot bead, 
which is the hit bead. The remaining EC displace¬ 
ment Xmax is decreased by x. 

4. Iterate by going back to step 2 until Xmax = 0. 

Then the next EC move is started. The relevant com¬ 
putational step in the EC moves is the evaluation of the 
admissible displacement. In a system consisting of un¬ 
bonded hard spheres, this is the distance to the bead 
hit first by the pivot bead while moving in the chosen 
direction. 

The EC algorithm can be adapted to spheres with p air- 
wise position-dependent interaction potentialsFor 
each move of the pivot bead in an EC chain, an en¬ 
ergy difference /S.E > 0 is drawn according to the Boltz¬ 
mann distribution. A displacement of the pivot bead 
that reduces the interaction energy is accepted (as in the 
standard Metropolis algorithm). A displacement increas¬ 
ing the energy is only partly executed, up to the point 
where the energy difference /S.E that has been drawn is 
reached or until the remaining EC displacement has been 
exhausted. 


Now we consider the general situation that the pivot 
bead has several pairwise interaction energies. For each 
interaction partner the energy difference /S.Ei then de¬ 
fines a maximal displacement of the pivot bead Xi < 
^max(AF^i). The largest possible displacement x of the 
pivot bead is the minimum of all which shall be real¬ 
ized for an interaction partner i.e., x = Xj = min^x^. 
The EC is then continued at bead j as next pivot bead. 
For hard sphere interactions, this algorithm reduces to 
the standard EC collision rule. 


Fig. I^b) shows an example for hard sphere polymers 
bonded by springs. The attempted EC displacements of 
the green bead are in one of two classes: (i) the beads do 
not collide along the path (for an EC move in direction 
vi) or (ii) the beads do collide along the path (for an EC 
move in direction V 2 ). In both cases, the energy stored 
in the bond reduces on the blue part of the trajectory, 
which, therefore, is always admissible, and increases on 
the red part, where the maximal admissible displacement 
is set by the “consumable” energy AE drawn from the 
Boltzmann distribution. Thus, the bond energy is only 
relevant for the maximal displacement if beads do not 
collide, because the other case is dominated by the hard 
sphere constraint. After displacing the green bead, the 
red bead becomes the pivot bead in both cases (i) and 
(ii). 


We prefer to choose the direction v of ECs randomly, 
which satisfies detailed balance. This can be relaxed, 
in principle, to other choices as discussed in Ref. [22] for 
hard sphere systems such that global balance is still sat¬ 
isfied. One particular simple choice, which can also be 
applied to the hard sphere polymers, is to start ECs only 
into three positive cartesian directions, which can gain a 
factor of approximately 2 in simulation speecP^ (essen¬ 
tially by simplifications in the collision detection). This 
simplification is not efficient, however, if combined with 
the parallelization scheme discussed in Sec. II A[ which 
decomposes the system into simulation cells and reflects 
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the ECs on simulation cell boundaries rather than reject¬ 
ing the whole EC move. Eor a cell decomposition with 
rectangular boundaries along cartesian directions, as it is 
usually used, EC moves started into cartesian directions 
will always reflect on themselves. 

A. Parallelization 

We use a parallelized version of this event-chain al¬ 
gorithm and refer to our earlier work for details of the 
parallelization.!^ As discussed there, the parallelization 
requires a spatial decomposition of the system (which is 
changed in every sweep to ensure ergodicity) into simu¬ 
lation cells. This limits the displacement of each sphere 
to its respective simulation cell. Eor non-bonded hard 
spheres this can be treated by reflection of non-admissible 
ECs at the cell boundaries. If a spring-triggered event oc¬ 
curs, where the bonded bead, which caused the event and 
would be the next pivot bead, is lying outside the current 
simulation cell, we proceed in a very similar manner, i.e., 
by reflection at the plane normal to the bond of the two 
participating spheres as illustrated in Eig.j^c): The gray 
bead is rendered immotile because of the currently cho¬ 
sen spatial decomposition into parallel simulations cells. 
Therefore, the EC cannot be transferred to the gray bead 
at the occurrence of a spring “collision”. Then, the pivot 
bead (green) does not change, and the propagation di¬ 
rection is reflected as if there was a wall normal to the 
bond. 

In this work we use a spatial decourposition scheme dif¬ 
ferent from a checkerboard partitiorP^ we use a rectan¬ 
gular tile-joint partition, where large tiles are separated 
by small joints (areas which contain spheres that cannot 
move). As discussed in Ref. |25l larger cells will lead to a 
more effective parallelization. 

B. Initial configurations 

The equilibration of polymer melts in simulations can 
be improved by generating initial configurations that are 
already representative of equilibrium configurations.!^ 
Erequently used strategies consi st in a slow compression 
of an equilibrated dilute solutiorP^ or a “push-off” pro¬ 
cedure, where the strongly repulsive steric interaction is 
switched on only after generating equilibrated configu¬ 
rations with a soft repulsive potential.!^ Eor the hard 
sphere polymer melt we propose an EC-based algorithm, 
which is conceptually similar to the slow push-off proce¬ 
dure in Ref. [21] for a Lennard-Jones melt. 

The flexible polymers in the equilibrated melt are ideal 
but acquire an effective stiffness. The effective stiffness is 
characterized by a finite value of (cos 0) {0 being the bond 
angle) .!^ Eor a long ideal chain of bond length cr, this re¬ 
sults in a mean-square end-to-end distance {R^){N) = 
cNa^ with a parameter c = (1 + (cos^))/(l — (cos^)) > 
IpH which depends on the short-range interaction be- 



FIG. 3. Mean-square internal distances R^(n)/ncr^ between 
two monomers with a chemical distance n along the chain for 
different initial condition generators and a long run simula¬ 
tion (thick black line). Small points are mean-square internal 
distances directly after setting up the phantom chains, cor¬ 
responding large points after increasing the bead size to cr, 
either by a “slow push-off” (red and blue curves) or by a 
“fast push-off” (green curve). The simulation parameters are 
N — 121, M = 500, L = 40cr (packing fraction rj = 0.495). 
We used ECs with a total displacement length i = 2a for the 
push-offs. 


tween polymer beads. Erom long-run simulation data, 
we find c ~ 1.9 for a melt of long hard sphere chains. 

In order to capture the effective stiffness already at 
the level of the initial configurations, we set up a system 
with randomly placed phantom polymers with vanishing 
hard sphere diameter and bond length <j, which we grow 
as non-reversal random walks by restric ting subsequent 
(unit) tangents to < cos(6>max)- Foi* an oth¬ 

erwise uniform distribution of bond vectors, this leads 
to (cos^) = COS^(^max/2). Wc choOSC ^max SUCh that 
{R^)/N 1.90(7^ holds in accordance with our long-run 

simulation data, see green and blue lines with small sym¬ 
bols in comparison to black line in Eig. 

We then introduce a finite excluded volume, but with a 
hard sphere diameter that is only a fraction of the target 
diameter a. This generates some “conflicts”, i.e., over¬ 
lapping spheres. We remove these conflicts by repeat¬ 
edly starting EC moves into different directions from the 
overlapping spheres only, until the conflicting overlap has 
been removed. In these ECs, we ignore pre-existent over¬ 
laps so that the EC will only be transferred to a bead the 
current pivot bead is not overlapping with. This proce¬ 
dure corresponds to locally “rattling” in the hard sphere 
system until enough space has been created around the 
overlapping bead to insert it. Once all conflicts for a 
given diameter are solved, we increase the diameter and 
continue iteratively until the target diameter a is reached. 
The iterative growth of sphere diameters (which we call 
“slow push-off” due to conceptual similarity with Ref. 
Ej) leads to a smaller change in the initial distribution 
of mean-square internal distances R^{n) between two 
monomers with a chemical distance n along the chain 
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FIG. 4. From left to right, we show EC moves with/without 
bead swapping. The currently selected pivot bead is high¬ 
lighted by a red halo. The EC direction is along the vector 
connecting the middle beads of the two “polymers”. The ini¬ 
tial EC (first column) is restricted by the hard sphere inter¬ 
action of the middle beads, thus we use the Metropolis algo¬ 
rithm when the two beads touch (second column), to evaluate 
whether to swap beads (third column, upper row) or transfer 
the EC (third column, lower row). Einally, the remainder of 
the displacement is performed (fourth column). 


(averaged over all chains), see curves with large symbols 
in comparison to corresponding curves with small sym¬ 
bols in Fig. For comparison, we also generate initial 
configurations by a fast increase of a (which we call “fast 
push-off” as in Ref. ETJl, see green curves in Fig.|3] 

Configurations after the push-off should exhibit in¬ 
ternal distances R^{n) as close as possible to the equi¬ 
librium result as found by a long simulation run, see 
black line in Fig. The initial configurations generated 
with slow push-off and the optimal value c ~ 1.9 (blue 
lines in Fig. are indeed similar to the long-run sim¬ 
ulation results. The fast push-off configurations (green 
lines) deviate with a maximum at intermediate N^ which 
takes a Iom time to equilibrate due to the slow reptation 
dynamics.^ Initial configurations generated from ideally 
flexible phantom chains (red lines) differ considerably. 

Using the slow push-off we can initialize systems at 
(in principle) any geometrically possible density with¬ 
out resorting to configurations that are far from equilib¬ 
rium (e.g., placing the beads on a lattice) in a reasonable 
amount of time (a couple of minutes wall tim^^ for the 
system parameters below). Even for a very dense sys¬ 
tem with T] = 0.63995, N = 120 and M = 275, we can 
generate initial configurations in 0{100h) wall time. For 
such dense systems, however, these are only valid config¬ 
urations, which are far from equilibrium because bonds 
are very elongated, and a thorough equilibration is still 
necessary. 


beads participating in the EC move. This results in a 
small collective translation of all beads participating in 
an EC cluster move with only small changes to the topol¬ 
ogy of entanglements. 

Topology changing MC moves such as the double¬ 
bridging mov^^ can spe ed up equilibration in polymer 
melts significantlyj ^ * ^^ * ^ ^ Here, we improve sampling with 
EC moves further by introducing an additional swap 
MC move, which can locally change the topology of en¬ 
tanglements. In contrast to the double-bridging move, 
which changes bonds, the swap move changes topology 
by changing bead positions. For this purpose, we modify 
the EC move so that the EC does not directly transfer to 
the next bead upon hard sphere contact but, instead, a 
swap of the two touching spheres is proposed, see Fig.[^ 
Such an additional swap move allows for a local change 
of entanglements. 

The EC swap move is accepted according to the stan¬ 
dard Metropolis algorithm. If the swap is rejected, the 
EC is transferred and the standard EC algorithm as de¬ 
scribed above is recovered. If it is accepted, the two beads 
are exchanged, and the EC continues with the same pivot 
bead. The example of a swap move in Fig. shows a sit¬ 
uation where it might be energetically favorable to swap 
beads. Note that in the absence of bonds all beads be¬ 
come indistinguishable, and the EC algorithms with and 
without swapping are identical up to book-keeping dif¬ 
ferences. 

The swap move is EC-specific: The EC automatically 
selects colliding pairs of beads for swapping; if the swap 
move is rejected, the EC move can continue without re¬ 
jection of the entire EC move. Moreover, detailed bal¬ 
ance is satisfied, and bead swapping can be included with 
very little computational overhead into the EC scheme. 
An analogous swap move in a standard MC algorithm 
needs to select pairs of beads such that the swap move 
has a reasonable acceptance rate (the particles have to 
be reasonably close). Moreover, the selection rule has to 
satisfy detailed balance (for example, simply proposing 
the nearest neighbor for swapping will lead to a violation 
of detailed balance). Therefore, there is no straightfor¬ 
ward analogue of the EC swap move in a standard MC 
simulation with local moves. 

Since the swap move locally changes topology and de- 
entangles polymers, the dynamics is no longer realistic if 
swap moves are applied. In particular, reptation dynam¬ 
ics is suppresses by swap moves (see numerical results 
below). On the other hand, this is the reason why swap 
moves can accelerate equilibration of the melt. 


C. Additional Bead Swapping 

Typical conformations in a dense melt consist of highly 
entangled polymers. In the dense limit the ECs become 
very long, i.e., the rather large displacement of an EC is 
distributed on a lot of very small displacements of many 


III. VALIDATION 

In order to verify our algorithm, we address structural 
equilibrium properties of chains in a polymer melt by 
investigating their typical shape as characterized by the 
moment of inertia tensoi!^ and the distribution of end- 
to-end distance.l^ These structural equilibrium quanti- 
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ties provide a detailed comparison across polymer melt 
simulation algorithms. The results of a comparison be¬ 
tween different MC and MD simulation algorithms are 
shown in the Appendix. We find quantitative agreement 
between the EC algorithm and standard MC and MD 
algorithms, and agreement with previous MC simulation 
results and theoretical predictions, where available. 

IV. PERFORMANCE AND DYNAMICS 


brated with short run simulations on the same worksta¬ 
tion with four CPUs for comparable results. 

We choose three different systems to investigate the 
influence of the occupied volume fraction r] and chain 
length N on algorithm performance (we use the same 
systems for the validation of equilibrium properties in 
the Appendix): 

1. System I: M = 400, N = 120 and rj = 0.390; 

2. System II: M = 500, N = 120 and r] = 0.490; 


For the comparison of the performance of different 
algorithms, we distinguish algorithms by whether they 
use the EC or standard Metropolis algorithm for (i) the 
hard sphere interactions and/or (ii) the bond spring in¬ 
teractions (“EC” for event chain, “MC” for standard 
Metropolis) and (iii) if the algorithm is executed par- 
allelly (par) or serially (ser) and (iv) if the swap move is 
used (swap). Accordingly, we introduce a naming scheme 
for algorithms where, for instance, “EC-MC-par” refers 
to a parallelized simulation, where hard sphere interac¬ 
tions are handled by the EC, springs handled by standard 
Metropolis algorithm, and the swap move is not used. 

We compare five different algorithms, namely EC-EC- 
par-swap, EC-EC-par, EC-EC-ser, EC-MC-par, and MC- 
MC-ser. This allows us to analyze the parallelization per¬ 
formance gains by comparing EC-EC-par/ser and check 
if we achieve the theoretical speed-up factor given by 
the number of processor cores. We do not parallelize 
the standard MC algorithm, because it was shown previ¬ 
ously that strong scaling is achievable.l^The comparison 
of MC-MC-ser/EC-EC-ser gives the algorithmic speed¬ 
up by using the event chain algorithm. The comparison 
EC-EC/MC-par demonstrates the advantage of using the 
event chain on the pair potential, i.e., the bonds. 

Additionally, we compare our results to those from 
MD simulations performed using the highly optimized 
LAMMPS package.l^ As hard spheres cannot be used 
in a force-based MD simulation, we compare with beads 
that are interacting via the repulsive part of standard 
Lennard-Jones potentials, whereas the bonds remain 
Hookean springs. The identification of the effective 
hard sphere radius of such soft Lennard-Jones spheres 
has been subject of prior research.l^ Our results show 
that identifying the Lennard-Jones length scale ctlj (de¬ 
fined by the zero of the full Lennard-Jones potential, 
^j(<^Lj) = Ulj(oo) = 0) with the hard sphere diame¬ 
ter a suffices for the purposes of this work. This comes 
with the advantage that we can use the same initial con¬ 
figurations (generated by our EC-based procedure) for 
the MD and MC evolutions. 


A. Diffusional dynamics and algorithm speed 

We compare the speed of different algorithms in terms 
of wall time. Since the simulations ran on different Cen¬ 
tral Processing Units (CPUs), all wall times were cali¬ 


3. System III: M = 250, N = 240 and r] = 0.490. 

This means the volume fraction 77 increases from System 
I to System II, whereas the polymer length N increases 
when going from System II to System III. 

In the following, we will compare the performance of 
these algorithms by the inter- and intrapolymer diffu- 
sional behavior of polymer chains using time-dependent 
MSDs of polymer beads. For a chain with bead positions 
Ti (i = I,..., A^) and center of mass R = ^ Ylf=i D, we 
measure the MSD functionsj ^^ l ^^ * 

gi{At) = ([rjv/ 2 (i +At) -rjv/ 2 (i)]^)t, (3) 

52(At) = ([(rjv/2(t + At) - R(t + At)) 

- {r^N/2{t)-'R{t))f)f (4) 


gi describes the diffusion of the middle bead including 
contributions from inter- and intrapolymer diffusion and 
g 2 the intrapolymer diffusion of the middle bead relative 
to the center of mass of the polymer. For both quantities, 
the average (.. .)t is an ensemble average and an average 
over time. 

In a polymer melt, the time evolution is governed by a 
sequence of crossovers 


gi{t) - < 


for t < Te 
for Te < t < Tr 
U/2 for TR<t <Td 

t for Td <t 


( 5 ) 


with three different crossover time scales: the entangle¬ 
ment time scale Te, the Rouse time scale tr^ and the dis¬ 
entanglement time scale r^^.^For all times scale t > Te^ 
reptation slows down the diffusional dynamics. The rel¬ 
ative MSD g 2 exhibits the same regimes as gi but is in¬ 
sensitive to center of mass diffusion. For t it ap¬ 

proaches a plateau value given by the radius of gyration 

i?2(7V) = iV-iEto((r-R)")- 

Any simulation dynamics achieving equilibration of in¬ 
trapolymer modes, will reach the plateau in the relative 
MSD ^ 2(^)5 beyond which intrapolymer fluctuations are 
equilibrated. We use the relaxation time to reach the 
plateau as a measure of simulation speed because it char¬ 
acterizes the equilibration performance of an algorithm 
on the scale of whole polymer chains. If the algorithm 
correctly describes the polymer melt dynamics on long 
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time difference At[s] 



FIG. 5. (a) Relative mean-square displacement ^ 2 (At) (log-log plot) for different algorithms against wall time (in s). g 2 

approaches a plateau value, which defines the polymer relaxation time scale, which serves as a measure for algorithm speed, 
(b) The mean-square displacements ^i(At) and ^2 (At) (log-log plots) for different algorithms with time rescaled to collapse all 
curves. Collapse is achieved for all algorithms except the one using swap moves. All simulations were performed for System I. 


time scales and exhibits Rouse, reptation, and chain dif¬ 
fusion dynamics as in Eq. (§, this relaxation time will 
coincide with the polymer disentanglement time Td. In 
Fig. E a), we show the wall time evolution of g 2 {t) for dif¬ 
ferent algorithms. Both MD (LAMMPS) and local MC 
dynamics f ollow Rouse dynamics with a t^/^-behavior for 
short times.f^^I^ Remarkably, we find such Rouse dynam¬ 
ics also for the cluster EC algorithm, even in the presence 
of of swap moves. All algorithms approach a plateau in 
the relative MSB g 2 {t). 

This also allows us to easily compare the performance 
of the algorithms and to determine a speed-up factor for 
each algorithm by rescaling time, i.e., shifting the double 
logarithmic curves such that the curves ^2 {t) coincide for 
long time scales close to the plateau. As a result, the sin¬ 
gle polymer relaxation time, which is identical to the dis¬ 
entanglement times Td if the algorithm exhibits all char¬ 
acteristic regimes of polymer melt dynamics, should be 
identical after rescaling. The resulting speed-up factors 
with respect to the standard local Metropolis algorithm 
MC-MC-ser are shown in Table HI 

If we use these speed-up factors for a linear rescaling 
of the time, the data for both MSB functions g 2 {t) and 
gi{t) and from all algorithms collapse onto two “master 
curves” as shown in Eig. [^b). The exception is the EC 
algorithm employing topology-changing swap moves. Be¬ 
cause this collapse includes the MB algorithm, this pro¬ 
vides evidence that both local MC dynamics and the clus¬ 
ter EC dynamics (in the absence of swap moves) evolve 
the system in a way that allows for an identification of 
“Monte Carlo time” (i.e., number of moves) with physical 
time. 

Eor the comparison in Table [l| we did not explicitly 
optimize the free simulation parameters like the total dis¬ 
placement length £ for an EC or the number of started 
ECs per sweep in a parallelized simulation (see Ref. |25| 
for a detailed discussion). Nevertheless, it is obvious that 
all algorithms clearly outperform the standard MC algo¬ 
rithm. Without parallelization, the EC-EC-ser algorithm 


algorithm 

System I 

System II 

System III 

MC-MC-ser 

1 

1 

1 

EC-EC-ser 

9 

7 

7 

EC-MC-par 

14 

17 

17 

EC-EC-par 

31 

25 

28 

EC-EC-swap-par 

230 

115 

423 

LAMMPS (par) 

330 

625 

770 


TABLE I. Comparison of relative speed-up compared to the 
standard local Metropolis algorithm MC-MC-ser for different 
system parameters and algorithms. Both parallel EC and 
LAMMPS simulations were performed using four cores. 


achieves speed-up factors up to 10 compared to the stan¬ 
dard MC algorithm (MC-MC-ser). The parallelization 
gives an additional speed-up factor of 3.5... 3.9 close to 
the theoretical limit of 4 given by the number of cores we 
used for the parallel simulation. We note that also stan¬ 
dard MC algorithms could be parallelized such that this 
additional parallelization speed-up factor is not specific 
to the EC algorithms. 

Bespite these speed-up factors for the EC algorithm, 
the LAMMPS MB simulation is still the fastest algo¬ 
rithm. Eor the comparison in in Table [T| we used a 
parallel version of LAMMPS running on four cores. We 
note that LAMMPS is under development since the mid 
1990^^ whereas our EC algorithm implementation, while 
adhering to general good practice rules for scientific com¬ 
putation, should still have room for optimization. In 
view of these preliminaries, the performance difference 
between the MB LAMMPS simulation and our fastest 
EC variant including swap moves seems very promising. 

Table [T] also shows that the EC-MC algorithm gains 
some efficiency with respect to EC-EC algorithms with 
increasing g. In such dense systems, the springs are com¬ 
pressed to a value close to their rest length a. Therefore, 
the rejection rate caused by the spring energy is rather 
low, such that the gain from the additional computa- 


















FIG. 6. Mean-square displacement ^i(At) (log-log plot) for 
different algorithms with time rescaled to collapse all curves 
for M = 20 polymers with length N = 500 (values for the 
spring constant k in units of kBTla‘^). 

tional effort in the rejection-free treatment of springs is 
small in denser systems. Since the disentanglement time 
Td ^ is strongly influenced by the chain length 
the efficiency of the swap algorithm increases with longer 
chains. MD performance does not decrease with density, 
whereas local MC and also EC performance depends on 
the displacement length which decreases with density. 
This explains the performance differences if the density 
T] is increased. 

The speed-up factors in Table [T| characterize algorithm 
equilibration times based on the polymer disentangle¬ 
ment time Td^ Alternatively, the autocorrelation time 
of the end-to-end vector c an been used to characterize 
equilibration times.^ ^^ * ^^ * According to Ref. [SU these 
equilibration time scales are comparable for moves not 
changing the chain topology; the disentanglement time 
Td is to be preferred if topology-changing moves are em¬ 
ployed that cannot relax density fluctuations (e.g., double 
bridging moves or our swap move). 

B. Reptation dynamics 

The good collapse onto master curves in Fig. [^b) sug¬ 
gests that we can observe the same regimes of polymer 
melt dynamics in the EC simulation as in a MD simula¬ 
tion, at least if no swap moves are employed. Therefore, 
we investigated whether also a regime of reptation dy¬ 
namics is observable with the EC algorithm. 

The reptation regime for Tq < t < tr rather hard to 
observe in simulations of shorter chains, and one expects 
a slightly increased exponent gi{t) ^ with 0.25 < x < 
0 dpi Yqy the chain lengths N = 120 used in Fig. [^b), 
the intermediate reptation regime is not clearly visible. 
In lattice MC simulations, evidence for an intermediate 
reptation-like regime with a considerably slower increase 
than has only been found in melts with long chains 
of length of N = 512.C^ 

Therefore, we also simulated a smaller system with less 


(M = 20) but longer polymers (A^ = 500) at = 0.298 
and measured gi{t) for EC-EC algorithms in compar¬ 
ison with MD simulations (using LAMMPS), see Fig. 

For these longer chains, the MD simulations show 
a much more pronounced intermediate regime of slowed 
down dynamics. We find that this regime is increasing for 
stiffer polymer springs {k = lOOksT /as compared to 
k = 30kBT/a‘^ in Fig. H): with a small probability, chains 
can still cross via thermally activated bond stretching, 
which becomes less probable for stiffer springs. For stiff 
springs {k = lOO/c^T/cr^), the MD simulations exhibit 
a reptation regime with a time-dependence gi{t) oc 
close to and in accordance with theoretical predic¬ 
tions and simulations in Ref. [32l 

The parallelized EC-EC simulations show exactly the 
same dynamical regimes in the MSD function gi{t) as the 
MD simulations, see Fig. In each EC move, all beads 
that would collide successively during a short time inter¬ 
val in a MD simulation are displaced at once. This gives 
rise to a MC dynamics which is effectively very similar to 
the realistic MD dynamics. For the EC-EC simulations, 
we considered a spring constant k = SOksT/a^ and, in¬ 
stead of stiff springs, hard sphere polymers bonded by 
tethers of maximal length = 1.4cr. In the tethered 
system, thermally activated bond crossing is absent sim¬ 
ilarly to a system with very stiff springs. Very similar 
to the stiff spring MD simulation, the tethered EC-EC 
simulation shows evidence of an intermediate reptation 
regime with gi{t) oc To our knowledge, this is the 
first off-lattice MC simulation, where clear indications of 
reptation dynamics could be observed. 

Fig.lH also shows, that the reptation regime is ab¬ 
sent as soon as we employ additional disentangling swap 
moves in accordance with our expectation. Swap moves 
can thus be used to accelerate equilibration by effectively 
“switching off” the slow reptation dynamics. 


V. CONCLUSION 

We introduced novel efficient off-lattice MC algorithms 
for the simulation of dense polymer melts of hard sphere 
polymers, which are based on event chain cluster moves 
previously known for hard sphere systems, see Fig. 
These EC cluster moves allow for a rejection-free treat¬ 
ment of the excluded volume interaction in the polymer 
melt. We generalize the algorithm to also handle the 
spring interactions in polymer bonds rejection-free. 

In addition, we introduce an efficient procedure to gen¬ 
erate initial configurations, which are representative of 
typical equilibrated configurations in polymer melts. Us¬ 
ing EC “rattling”, we can generate initial configurations 
up to very high packing fractions (up to = 0.63995). 

We parallelize the event chain Monte Carlo algorithm 
and suggest additional local topology-changing swap 
moves, see Fig. to further increase simulation speeds 
in melts. 

We validated the EC algorithm by comparing equi- 











9 


librium structural properties. In the Appendix, we show 
results for the polymer shape (Fig.[^ and the end-to-end 
distance distribution (Fig. |^, which are in quantitative 
agreement with other MD and MC simulation techniques. 

We assessed the performance of the EC algorithm by 
measuring its equilibration speed using the relative MSD 
function g 2 (t) of a polymer bead in the middle of a poly¬ 
mer with respect to the polymer center of mass, see Fig. 

This allows us to define a polymer relaxation time, 
which is specific to the algorithm and represents a mea¬ 
sure for its equilibration speed. We find that the combi¬ 
nation of EC moves and parallelization can increase MC 
simulation speeds by factors up to 30. If also swap moves 
are employed, MC simulation speeds become comparable 
to optimized MD simulations that we performed with the 
LAMMPS package for comparison. 

Without swap moves, the dynamics of the EC algo¬ 
rithm appears to be very similar to MD dynamics. A 
simple rescaling of simulation times can collapse MD and 
EC simulation dynamics, see Eigs. [^and[^ The collective 
dynamics generated by the EC moves, which essentially 
displace all beads coherently that collide successively in 
a short time interval in a MD simulation, appears to be 
very similar to the MD collision dynamics. 

Accordingly, in the absence of swap moves, the EC 
algorithm exhibits all dynamical regimes expected for 
polymer melts, i.e.. Rouse, reptation, and chain diffusion 
dynamics. In particular, we can identify an intermedi¬ 
ate reptation regime with a MSD function gi{t) oc 
close to in simulations of a system with long chains 
(A^ = 500), see Eig. To our knowledge, this is the 
first off-lattice MC simulation, where reptation dynam¬ 
ics could be observed. If topology-changing swap moves 
are used, which disentangle polymer chain, reptation dy¬ 
namics is absent in the EC algorithms. 

Although we only presented results for the most sim¬ 
ple case of a melt of flexible polymers with no interpoly¬ 
mer interaction other than excluded volume, the added 
value of EC algorithms should persist in more complex 
systems. Eor (bond) interactions that are not pair inter¬ 
actions, e.g., bending energies, reject ion-free sampling in 
the way presented here is not possible. We have already 
shown in a previous work,l^ however, that such bending 
energies can still be treated by proposing moves that are 
compliant with the hard sphere constraint by using ECs 
and than accepting (or declining) this move according to 
the standard Metropolis algorithm. 
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Appendix A: Validation 


1. Moment of inertia tensor 

The shape of a polymer in a dense melt is ellipsoidal 
rather than spherical, which can be shown by the distri- 
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EIG. 8. Distribution of the end-to-end distance Wl{R) and 
in comparison with eq. (A4) for ideal chains. Eor clarity the 
results for System I are shifted by an offset. Eor short dis¬ 
tances Wl{R) exhibits oscillations. The inset shows the ratio 
Wl{R)/ g{r) (red points) with the pair correlation g{r), which 
shows no oscillations. 



EIG. 7. (a) Bead distribution along the three principal axes 

of the moment of inertia tensor (for System I). (b) Snapshots 
of polymer configurations. The left column shows a single 
configuration, the right column an overlay of 50 configurations 
to visualize the distribution of the beads. Each row shows the 
projection to a plane spanned by two distinct principal axes. 
The ellipsoids have the same moment of inertia tensor as the 
polymer. 


bution of beads with respect to the center of mass in the 
coordinate system which is given by the eigenvectors 
of the moment of inertia tensor, 


(I)ii = - rk,irk,j), (Al) 

k 


of a polymer .1^ The sum runs over all beads of a polymer, 
where rk,i denotes the ith component of the kth bead 
coordinate. 

Following Ref. 123 we can use the eigenvalues Ii < 
I 2 < h of the moment of inertia tensor to characterize 
the shape in terms of its asphericity b = |(/i + 12 ) — /s, 
acylindricity c = /i — / 2 , and shape anisotropy = 
4(1 — 3 ( 12/3 + /3/1 + /i/ 2 )/(Tr I)^. Additionally, there 
exist several analytical predictions for an infinite freely 
jointed chain,!^ 


and foil^ 


lim 

N^oo 


lim 

N^oo 


(J1+J2-J3) 

(TrI) 

(J1-J2 + J3) 

(IVI) 


lim 

N^oo 


{-h + /2 + / 3 ) 

(TrI) 


= 0.754, 


= 0.175, 


= 0.0646, 


(A3) 


which can be tested. 

In Fig. we compare the distribution of beads of one 
polymer in the system spanned by the eigenvectors of 
the moment of inertia tensor for all algorithms (for Sys¬ 
tem I). The different widths of the distributions along 
the three principal axes of the moment of inertia ten¬ 
sor in Fig. I^a) implies that polymers in the melt have 
an ellipsoidal shape. The distribution along the largest 
eigenvalue axis is bimodal corresponding to an additional 
dumbbell-like shape in this direction. The agreement be¬ 
tween all simulation algorithms is excellent. Our results 
also agree with MC simulation results in Ref. [28l In Fig. 
IZlb), we visualize the actual shapes of polymers demon¬ 
strating the prolate shape of a polymer. Snapshots in the 
first two rows confirm the the dumbbell-like shape with 
a minimum in the bead distribution along the largest 
eigenvalue axis. 

In Table|nl we compare the shape descriptors from our 
Monte-Carlo schemes and LAMMPS with the theoretical 
expectations ( |A2[ ) and (A3). All results coincide very 
well. 


2. Distribution of the End-to-End Distance 


(462+3c2) 2 

(Trl)^ = 3 


The distribution of the end-to-end distance W{R) for 
(A2) ideal chain with a mean-s quare end-to-end distance 

{R‘^){N) = cNa‘^ (see section IIB for the definition of 
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Sys.I 

MC-MC-ser 

Sys.II 

Sys.III 

EC-EC-swap-par 

Sys.I Sys.II Sys.III 

Sys.I 

LAMMPS 

Sys.II 

Sys.III 

theo. 

(K^> 

0.390 

0.388 

0.387 

0.399 

0.397 

0.395 

0.402 

0.397 

0.395 

0.42 

(/i + / 2 -/ 3 )/(Tr(I)) 

0.763 

0.758 

0.757 

0.765 

0.764 

0.766 

0.765 

0.767 

0.763 

0.754 

{h -I 2 + l 3 }/{Tr{l)} 

0.173 

0.178 

0.177 

0.172 

0.173 

0.172 

0.172 

0.170 

0.173 

0.175 

{-h +I 2 + / 3 )/(Tr(I)> 

0.0642 

0.06473 

0.06540 

0.06284 

0.06310 

0.06251 

0.06225 

0.06302 

0.06427 

0.0646 

(462 + 3c2)/(Tr(I))2 

0.653 

0.621 

0.629 

0.640 

0.638 

0.672 

0.646 

0.643 

0.633 

0.667 


TABLE IL Comparison of shape descriptors from MC algorithms, LAMMPS, and theoretical expectations. For clarity, we only 
show the results from MC-MC-ser, EC-EC-swap-par and LAMMPS algorithms; the other variants do not differ significantly. 


the stiffness parameter c) is approximately given by a 
Gaussian distributiorP^H 

W(RHR = exp (- 

(A4) 


Therefore, the ideality of chains in a polymer melt can be 
tested by comparing simulation results for the distribu¬ 
tion of the end-to-end distance with the Gaussian expec¬ 
tation (A4), see Fig.[^ The agreement with the Gaussian 


expectation is indeed good, apart from an oscillating be¬ 
havior at small distances R. These oscillations can be 
explained by the influence of the pair correlation func¬ 
tion g{r) characterizing the additional local liquid-like 
ordering of neighboring polymer beads. These oscilla¬ 
tions are in quantitative agreement with g{r)WL{R)dR^ 
where we determined the pair correlation g{r) of beads 
in the polymer melt numerically. 

Also the agreement among the results for different sim¬ 
ulation algorithms in Fig. is very good. Only the stan¬ 
dard serial MC-MC algorithm shows deviations because 
of its long equilibration times. 










